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Abstract. A measurement of the cosmic-ray spectrum for energies exceeding 4xl0 18 eV 
is presented, which is based on the analysis of showers with zenith angles greater than 
60° detected with the Pierre Auger Observatory between 1 January 2004 and 31 Decem¬ 
ber 2013. The measured spectrum confirms a flux suppression at the highest energies. 
Above 5.3xl0 18 eV, the “ankle”, the flux can be described by a power law E ± 7 with in¬ 
dex 7 = 2.70 ±0.02 (stat) ±0.1 (sys) followed by a smooth suppression region. For the energy 
(E s ) at which the spectral flux has fallen to one-half of its extrapolated value in the absence 
of suppression, we find E s = (5.12 ± 0.25 (stat)^} ^ (sys))xl0 19 eV. 
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1 Introduction 

The Pierre Auger Observatory is to date the largest detector of air showers induced by the 

ultra-high energy cosmic rays (UHECRs). It is a hybrid detector combining an array of 

surface detectors (SD), that measures the particle densities in the shower at the ground, 

and a fluorescence detector (FD), that captures the ultraviolet light emitted by nitrogen 

as showers develop in the atmosphere. The SD [1, 2], spread over an area of 3000 km 2 , is 

composed of a baseline array, comprising more than 1600 water-Cherenkov detectors placed 

in a hexagonal grid with nearest neighbours separated by 1500 m. The FD [3] comprises 

24 fluorescence telescopes distributed in 4 perimeter buildings which view the atmosphere 

over the array. The site is located near Malargiie, Province of Mendoza, Argentina, at an 

altitude of about 1400 m above sea level and at an average geographic latitude of 35.2° S. The 

redundancy provided by the two detection techniques has proved to be extremely valuable for 

a wide range of applications and has improved the performance of the Observatory beyond 

expectations. 

The data gathered at the Pierre Auger Observatory with zenith angles less than 60° 

have provided a measurement of the UHECR spectrum [4] with unprecedented statistics. 

The technique developed exploits the large aperture of the SD, operating continuously, as 
well as the calorimetric measurement of the energy deposit obtained with the FD which is, by 
contrast, rather limited in duty cycle to clear nights without moonlight (13%). A parameter 
quantifying the shower size is obtained from the SD. This parameter is then calibrated using 
the energy inferred from the calorimetric FD measurement for a sub-sample of the events 

(hybrid events) which are detected and reconstructed simultaneously with both techniques [4] . 
This allows measurements of the energy spectrum with an energy estimation weakly reliant 
upon shower simulations. The spectrum obtained is consistent with the Greisen and Zatsepin- 
Kuz’min (GZK) suppression [5, 6]. The observation of this strong flux suppression was first 
reported by HiRes [7]. In addition, the spectrum has been independently measured using 
hybrid events which are detected with the fluorescence technique and at least one particle 
detector [8]. The latter measurement has also been combined with the SD spectrum to 
obtain a spectrum extending to lower energies [8, 9], which has confirmed the flattening 
of the spectral slope at about 5xl0 1 * * * * * * * * * * * * * * * * 18 eV, often referred to as the “ankle”. Similarly the 
Telescope Array (TA) experiment [10, 11], having an array of scintillators spread over an area 
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of 700 km 2 , together with fluorescence telescopes, has provided an independent measurement 
of the spectrum using showers with zenith angles less than 45° in the Northern Hemisphere, 
which clearly displays the GZK steepening [12]. Although the two spectra are in agreement 
between the ankle and 5xl0 19 eV, there are discrepancies in the suppression region which 
are being scrutinized by a joint working group of both Collaborations. 

The choice of deep-water Cherenkov detectors for the Pierre Auger Observatory was 
made to give substantial solid angle coverage for zenith angles exceeding 60° thus enabling 
detection of showers over a much broader angular and declination range than is possible in 
arrays using thinner detectors. The showers with large zenith angles traverse much larger 
atmospheric depths than than showers nearer the vertical and, as a result, are characterised 
by the dominance of secondary energetic muons at the ground, since the electromagnetic 
component of the shower is largely absorbed before reachig ground. Inclined showers have 
been observed since the 1960s and have stirred much interest [13], but the first reliable 
methods for event reconstruction are relatively recent [14, 15]. These events provide an 
independent measurement of the cosmic-ray spectrum. In addition, since inclined showers 
are mainly composed by muons, their analysis has been shown to be particularly useful 
in determining the muon content of showers for comparison with predictions from shower 
models [16] . Inclined showers are also used to extend the accessible fraction of the sky covered 
in analyses of the distribution of the arrival directions of cosmic rays detected at Auger [17, 18] 
and constitute the main background for searches for ultra-high energy neutrinos using air 
showers [19]. 

In this article we present the measurement of the cosmic-ray spectrum derived from 
events with zenith angles between 60° and 80° detected with the Pierre Auger Observatory 
in the time period from 1 January 2004 to 31 December 2013. 

2 Estimation of the shower energy 

Since inclined showers are dominated by muons that are strongly affected by the geomagnetic 
field, they require a specific reconstruction procedure that differs from that used for events 
with zenith angles less than 60°. The method to estimate the energy of an inclined event 
recorded by the SD array, described in full detail in [15], is summarised in the following. 

The arrival direction (9 , eft) of the shower is obtained by fitting the start times of the 
signals recorded at the triggered stations to a plane front corrected for small time delays 
associated with the muon propagation as modelled in [20]. Once the shower geometry is 
established, the intensities of the signals are used to quantify the shower size and the impact 
point of the shower axis on the ground. The method is based on the observation that the 
shape of the muon distribution is approximately universal for a given shower direction and 
that only the overall normalisation of the muon distribution depends on the shower energy 
and primary mass [15, 21]. It has also been shown that the lateral shape of the muon density 
profile is consistently reproduced by different hadronic models and software packages used for 
air-shower simulations. These universal characteristics allow one to model the muon number 
density as a function of the position at the ground (muon distribution) as: 

Pfi(r) = Nig p M ,ig(r; 0, <j>). (2.1) 

Here IVig is a measure of the shower size and is the relative normalisation of a particular 
event with respect to a reference muon distribution, p Ati ig(r; 9, (ft) [22, 23], and p^ig is conven¬ 
tionally chosen to be the average muon density for primary protons with an energy of 10 19 eV 
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simulated using QGSJetII- 03 [24] as the hadronic interaction model. Nig is thus expected 
to be correlated with the shower energy, and independent of the shower arrival direction. 
Therefore, Nig can be used as the energy estimator of inclined events recorded at the SD. 

Nig is estimated by fitting the measured signals at the SD stations to the expected 
muon distribution for a given arrival direction, using a maximum-likelihood method based 
on a probabilistic model of the detector response to muon hits obtained from Geant4 [25] 
simulations with the Off line framework [26] of the Pierre Auger Observatory. A residual 
electromagnetic component of the signal stemming from the muons themselves, mainly from 
decay in flight, is taken into account based on simulations (typically amounting to 20% 
of the muon signal) [27]. To ensure a good reconstruction, only events well-contained in 
the SD array are selected. This fiducial trigger requires that the reconstructed core of the 
shower falls within an active cell, defined as the elemental hexagon cell 1 surrounding an 
operational station when the six neighbouring stations in the regular hexagonal pattern are 
also operational, at the time of the event. 

Although the actual value for Nig depends upon the particular choice of composition and 
hadronic model made for the reference distribution, the absolute energy calibration is inferred 
from a sub-sample of hybrid events measured simultaneously with the FD and SD that are 
used to calibrate the shower size Nig with the calorimetric shower energy measured with the 
FD, Efd . This ensures that most of the uncertainties associated with the unknown primary 
composition in data and hadronic model assumed, as well as many uncertainties associated 
with the reconstruction, are absorbed in this robust and reliable calibration procedure based 
only on data. 

The FD provides a nearly calorimetric, model-independent energy measurement, be¬ 
cause the fluorescence light is produced in proportion to the electromagnetic energy deposited 
by the shower in the atmosphere. The longitudinal profile of the energy deposit, dE/dX, is 
determined from the signals at the triggered pixels of the telescope cameras after taking into 
account the separate fluorescence and Cherenkov light contributions, as well as light atten¬ 
uation and dispersion in the atmosphere [29] . The electromagnetic energy released by the 
shower in the atmosphere is obtained by fitting the longitudinal profile d E/dX to a Gaisser- 
Hillas (GH) function [30] and integrating over the X-range. The total primary energy, Epu, 
is then derived by adding the so-called “invisible energy” carried into the ground by high- 
energy muons and neutrinos, which is estimated using an unbiased and model-independent 
method [31]. The shower-energy estimated with the FD has a total systematic uncertainty 
of 14% [32], 

For the calibration, only events with zenith angles greater than 60° independently trig¬ 
gered and reconstructed by the FD and SD are considered. The FD measurements must pass 
quality cuts designed to select high-quality longitudinal profiles observed under good atmo¬ 
spheric conditions, including the condition that the depth of the shower maximum, X max , 
is within the field of view (FOV) of the telescopes. The latter condition together with the 
limited FOV of the FD may introduce a different selection efficiency for different primary 
masses. To avoid this effect, a strict cut on the slant depth range observed by the telescopes 
is also applied to ensure that the FOV is large enough to observe all plausible values of the 
depth of the shower maximum for the geometry of each individual shower. See [15] for a full 
description of the selection criteria. Finally, only events with energies above 4xl0 18 eV are 
accepted to ensure practically full trigger efficiency (see Section 3). 

1 The elemental hexagon cell is the region with vertices at the barycentre of each of the six triangles around 

the central station, with area given by D 2 \/ 3/2 where D = 1.5km is the array spacing [28]. 
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Figure 1. Correlation between the shower size parameter, Nig, and the reconstructed FD hybrid 
energy, Epu, for the selected hybrid data with 9 ^ 60° used in the fit. The uncertainties indicated 
by the error bars are described in the text. The solid line is the best fit of the power-law dependence 
-Esd = A ( Nig) B to the data. The corresponding ratio distribution of the SD energy Esb to the FD 
energy FIfd is shown in the inset. 


Here we present the calibration results using inclined events with zenith angles between 
60° and 80° recorded from 1 January 2004 to 30 September 2013, increasing the data sample 
by about 16% with respect to the one used in previous analyses [15, 16]. A total of 255 
hybrid events are selected. 

The correlation between the energy estimator Nig and the calorimetric hybrid energy 
.Efd is well described by a simple power-law function: 

N 19 = A'(E FD /10 19 eV) B '. (2.2) 

The fit is based on a tailored maximum-likelihood method [33] that includes both the uncer¬ 
tainties of IVig and £fd without relying on approximations. For each event the uncertainty 
in Nig corresponds to the reconstruction uncertainty as obtained from the fit to the reference 
muon distribution. Shower-to-shower fluctuations of Nig also contribute to the spread and 
are taken into account. The relative uncertainty associated with shower-to-shower fluctua¬ 
tions is assumed to be constant in the energy range of interest and fitted to the data as an 
extra parameter. The uncertainty in F?fd assigned to each event comes from the propagation 
of both the statistical uncertainties of the reconstruction procedure (associated with the re¬ 
construction of the shower geometry, the GH fit to the longitudinal profile and the correction 
for the invisible energy) and the uncertainties of the atmospheric parameters [32]. 

The resulting fit to the selected data is illustrated in figure 1. By inverting the fitted 
function, the SD energy estimate is given as Esb = A (N\g) B and the corresponding calibra¬ 
tion parameters are A = (5.701 ±0.086)x 10 18 eV and B = 1.006T0.018. Although statistical 
uncertainties of the calibration constants A and B affect the SD energy scale, their contribu¬ 
tion is small (from 1.5% at 10 19 eV to 4.5% at 10 2 °eV), decreasing as the number of events 
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increase. There is also an uncertainty of ~2% attributed to the different angular distributions 
of the hybrid events and the full inclined data set used to calculate the spectrum. The main 
contribution to the systematic uncertainty comes from the overall systematic uncertainty of 
14% from the FD energy scale. By adding the uncertainties described above in quadrature, 
the total systematic uncertainty in Egn ranges between 14% at 10 19 eV and 17% at 10 2 ° eV. 

The resolution in Ego is computed from the distribution of the ratio of A (Nig) B / Epp> 
for the hybrid events used for the calibration, assuming a fixed FD energy resolution of 
7.6% [32] (method based on [34]). The resulting average SD energy resolution is estimated 
to be (19 ± 1)% (see inset in figure 1). It is dominated by low-energy showers due to the 
limited number of events at the highest energies. 

3 Efficiency and exposure 

The final step in measuring the energy spectrum is a precise determination of the exposure 
for the observations. The exposure involves the time integral of the aperture, that is, the 
integral of the instantaneous effective area of the SD over solid angle and observation time. 
This integral is subsequently weighted by the trigger and event selection efficiency, which 
depends on the characteristics of the shower such as the nature of the primary particle, its 
energy and arrival direction. For low energies, the efficiency is smaller than one due to the 
trigger and the selection procedures (details in [15, 28]) and can in practice be a source of 
large uncertainty. It is thus advantageous to limit the spectral measurements to energies at 
which the array is fully efficient, that is, when the effective area of the SD coincides with the 
geometrical one. 

The inclined SD data set used to measure the cosmic-ray spectrum is composed of 
events selected to have zenith angles between 60° and 80° and falling inside an active region 
of the array to guarantee that no crucial part of the shower is missing. Hence, the trigger 
and event selection efficiency of the SD array for inclined showers is the probability that a 
cosmic-ray shower reaching an active region of the array induces a trigger, is selected and 
finally reconstructed. 

The energy for full detector efficiency is determined using hybrid events (detected with 
the FD and having at least one triggered SD station). To avoid biases from the primary 
composition, the same data selection criteria as for the energy calibration are applied (see 
Section 2). Additionally, it is required that the shower core reconstructed with the FD 
technique falls within an active cell of the SD array. Assuming that the detection probabilities 
of the SD and FD detectors are independent, the average detection efficiency as a function of 
calorimetric shower energy reconstructed with the FD can be estimated from the fraction of 
hybrid events that trigger the SD, are selected and finally reconstructed [15]. The resulting 
efficiency is shown in figure 2 as a function of both shower energy E and shower size parameter 
Nig. The conversion of IV 19 to energy has been made using the calibration procedure based 
on data discussed in the Section 2. The efficiency evaluated with the hybrid approach is 
subject to relatively large uncertainties, indicated as error bars, which are due largely to the 
limited number of events. The trigger efficiency of the SD array is found to be almost fully 
efficient (>98%) for energies above 4xl0 18 eV or for N±g > 0.7. This value will be used as 
our estimate of the energy for full efficiency. 

The energy for full efficiency obtained with the hybrid apprach is cross-checked using 
full shower and detector simulations. The sample used consists of about 20 000 proton show¬ 
ers simulated with CORSIKA [35] using QGSJetII-04 [36] with zenith angles isotropically 


- 5 - 




log l0 (E/eV) 


Figure 2. Trigger and event selection efficiency of the SD array for showers with zenith angles 
between 60° and 80° as a function of shower energy derived from the hybrid data (circles) and from 
the Monte Carlo simulated showers (squares). The error bars indicate the statistical uncertainty (the 
68 % probability contour). 


distributed between 60° and 80° and energies ranging from log 10 (£'/eV) = 18 to 20, in steps 
of 0.5 (with a spectral index 7 = 1 in each sub-interval). To generate the simulated events, 
these showers subsequently underwent a full simulation of the detector response, within the 
Off line framework, with random impact points in the SD array. Only showers with impact 
points in an active region of the array are considered for this analysis. For each event, the 
full trigger and event selection chain are applied. 

The trigger efficiency for inclined showers depends mainly on the total number of muons 
at the ground for a given arrival direction, which is proportional to the reconstructed shower 
size parameter N\g. A recent study [16] has reported that the measured muon number in 
inclined showers at ground level exceeds expectations obtained with simulations and various 
hadronic interaction models (even when assuming a pure iron composition). This implies 
that simulations and data have a different Nig—E correlation, that is, the MC energy is not 
directly comparable to Ep d for the same muon content. To carry out a comparison between 
using hybrid data and simulations, we correct for the muon deficit observed in simulations 
by converting the shower size of each simulated event into energy using the same calibration 
procedure applied for the data. The resulting efficiency as a function of both shower energy 
and shower size parameter is also shown in figure 2. The trigger and selection become almost 
fully efficient at energies above 4xl0 18 eV in agreement with hybrid data. It has been tested 
with simulations that above this energy, the detector is fully efficient independently of the 
particular choice of mass composition and hadronic model. 

The choice of a fiducial trigger based on active hexagons allows us to exploit the regu¬ 
larity of the array [28] . The geometrical aperture of the array is obtained as a multiple of the 
aperture of an elemental hexagon cell. Above the energy for full efficiency, the detection area 
per cell is 1.95 km 2 , which results in an aperture of a ce u = 1.35 km 2 sr for showers with zenith 
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angles between 60° and 80°. The calculation of the integrated exposure over a given period 
of time simply amounts to counting the operational cells as a function of time, N ce \\(t), and 
integrating the corresponding aperture, N ce n(t) x a ce ii, over time. The overall uncertainty on 
the integrated exposure is less than 3% [28]. 

In January 2004 the initial area spanned by the array was about 34 km 2 and it rose 
steadily until August 2008 when over 1600 surface detectors were in operation. The calcu¬ 
lation of the surface area is the same as that used for the spectrum measured with events 
at zenith angles less than 60° and has been done numerically, making use of the trigger rate 
at each station [28]. The integrated exposure for showers with zenith angles between 60° 
and 80° amounts to (10 890 ± 330)km 2 sryr in the time period from 1 January 2004 to 31 
December 2013 (which corresponds to 29% of the exposure for showers with 0<6O° in the 
same period). This calculation excludes periods during which the array was not sufficiently 
stable, which add up to a small fraction of the total time (of the order of 5%). 


4 Spectrum 

Inclined events recorded by the SD in the time period between 1 January 2004 and 31 De¬ 
cember 2013 were analysed and calibrated with the procedure outlined above. The resulting 
data set consists of 254 686 events that fell on the active part of the array with zenith angles 
in the interval between 60° and 80°. To avoid large uncertainties due to the trigger and event 
selection efficiencies, only events having energies greater than 4xl0 18 eV are considered for 
the spectrum calculation, reducing the sample to 15 614 events. 

The differential flux of cosmic rays at a certain energy E , J(E), is obtained by dividing 
the energy distribution of the cosmic rays by the accumulated exposure in the corresponding 
zenith-angle interval for the energies in the range of full trigger efficiency. Figure 3 shows this 
ratio for the selected events (also displayed as a function of the equivalent shower size IV 19 ). 
This is a raw distribution since effects of the finite resolution of the SD energy measurement 
have not yet been taken into account. Note that this figure also shows the ratio below 
4xl0 18 eV where it is not expected that the flux will be reproduced accurately because of 
the energy dependence of the triggering efficiency (see figure 2 ). 

A correction must be applied to account for the effect of the resolution in the energy 
determination, responsible for a bin-to-bin event migration. The number of events in a given 
bin is contaminated by movements of reconstructed energies from the adjacent bins. For 
an energy spectrum which is steeply falling, upward fluctuations into a given bin are not 
completely compensated by other fluctuations from the other direction, and the net effect is 
an overestimate of the flux when this correction is not applied. Due to this, the measured 
spectrum is shifted towards higher energies with respect to the true spectrum. To correct for 
this, a forward-folding approach is applied. Monte Carlo simulations are used to determine 
the resolution of the SD energy estimator based on the shower size parameter, IV 19 , for 
different assumptions on the primary mass and hadronic interaction models. Then the average 
resolution is converted to the SD energy resolution using the same energy scale as obtained 
from the data. In addition, intrinsic fluctuations (also called shower-to-shower fluctuations) 
contribute to the bin migration. In this case, these fluctuations are modelled with a normal 
distribution that has a constant relative standard deviation, a(N±g)/Nig, and are estimated 
in the data-calibration procedure as explained in [15], resulting in a(Nig)/N\g = (14 ± 1)%. 
From the SD energy resolution and intrinsic fluctuations, a bin-to-bin migration matrix is 


-7 - 


N , 9 



log i 0 (E/eV) 


Figure 3. Uncorrected energy spectrum of the cosmic rays derived from inclined events with zenith 
angles 60° ^ 6 ^ 80° in terms of energy and shower size. Only statistical uncertainties are shown. 
The number of events in each of the bins above 4xl0 18 eV (vertical dashed line) is given above the 
points. 


derived. The matrix is then used to find a flux parameterisation that matches the measured 
data when forward-folded. 

We assume a common parameterisation for the spectrum given a power law below the 
ankle and a power law plus a smoothly changing function above, 


J{E) oc E ~ 71 
J(E) oc E ~ 72 



j E <C -f^ankle; 


; E 5? E. 


ankle 9 


(4.1) 

(4.2) 


where E s is the energy at which the flux falls to one-half of the value of the power-law 
extrapolation of the intermediate region and Ay gives the increment of the spectral index 
beyond the suppression region - '. This parameterisation was convolved with the migration 
matrix and the resulting flux was fitted to the measured spectrum, using a binned maximum- 
likelihood approach assuming Poisson statistics. As Tinkle is very close to the saturation 
energy for inclined showers, neither i?ankie nor 7 i can be reliably established from the data. 
These parameters are relevant for the unfolding of the lower-energy part of the spectrum. 
They were fixed to the values obtained in the spectrum reported in [9], 71 = 3.23 ± 0.01 ± 

2 Note that Ay = l/ln(W c ) in the equivalent expression used in [8], where W c defines the width of the 
transition region at the suppression. 
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Figure 4. Correction factor applied to the measured spectrum to account for the detector effects as 
a function of the cosmic-ray energy. The uncertainty due to the energy resolution is shown with the 
dark band, and that due to shower-to-shower fluctuations with the blue band. The total uncertainty, 
also including the uncertainty of the fit, is shown by the light band. 


0.07 (sys) and © an ki e = (5.25 ± 0.12 ± 0.24 (sys))xl0 18 eV , while the other parameters were 
left free. The convolved flux is finally divided by the input flux to obtain the correction 
factors which are in turn applied to the measured binned spectrum. 

The correction factor resulting from the fit of the assumed parameterisation is shown in 
figure 4. The uncertainty from the fit was obtained by propagating the covariance matrix of 
the fitted parameters into the correction function. A systematic uncertainty is attributed to 
possible variations of the shower-to-shower fluctuations, which can arise if the mass compo¬ 
nent changes and because of the changing -A max over the energy range of interest. Although 
recent results [16] favour a transition from lighter to heavier elements in the energy range 
considered, here we assume a conservative scenario where the relative fluctuations are allowed 
to vary between 4% (corresponding to a pure iron composition) and 21% (corresponding to 
a pure proton composition) over the full energy range. Finally, the propagated uncertainty 
in the average resolution of the energy estimator N\g, estimated to be on the order of 10%, 
was also included. The uncertainty in the correction factor that arises from the uncertain¬ 
ties in the ankle parameters (.©ankle and 71 ), fixed to the values reported in [9], is negligible 
(<1% below E s ). In figure 4, in addition to the total uncertainty (light band), the last two 
uncertainty components are shown separately to illustrate that the main contribution is the 
systematic uncertainty due to variations of the shower-to-shower fluctuations. 

The spectrum, corrected for energy resolution, is shown in figure 5. The error bars 
represent statistical uncertainties only. The flux is multiplied by © 3 to better present its fea¬ 
tures, a flat region above 4xl0 18 eV up to the steepening at energies above about 4xl0 19 eV. 
The light shaded boxes indicate the total systematic uncertainties (less than ~40% up to 
4xl0 19 eV, and then increasing up to ~200% for the highest energy bin), which include 
the uncertainty in the calibration parameters propagated to the flux, a global uncertainty 
derived from the SD exposure calculation (3%), the uncertainty arising from the unfolding 
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Figure 5. Energy spectrum of the cosmic rays, corrected for energy resolution, derived from inclined 
data with zenith angles 60° ^ 9 ^ 80°. The error bars represent statistical uncertainties. The light 
shaded boxes indicate the total systematic uncertainties. For illustration purposes, the systematic 
uncertainties excluding the uncertainty from the energy scale are also shown as darker boxes. The 
effective number of events, after correcting the flux for energy resolution, is given above the points. 


process, and the global systematic uncertainty of the FD energy scale (14%) from the hybrid- 
calibration procedure. In figure 6 these separate contributions to the systematic uncertainties 
in the derived flux are shown as a function of the cosmic-ray energy. 

5 Discussion 

The characteristic features of the spectrum have been quantified by fitting the model given by 
eqs. (4.1) and (4.2), which assumes a spectrum with a sharp ankle and a smooth suppression 
at the highest energies, to the unfolded spectrum. The result of the best fit is shown as a 
solid line in figure 7. Another approach is to describe the spectrum with three power laws 
separated by two breaking points as illustrated by a dashed line in figure 7. The first model 
improves marginally the description of the abrupt change of slope observed at higher energies. 
The spectral parameters from the best fits to the data are given in table 1, quoting both 
statistical and systematic uncertainties. 

Note that, as explained above, the selected energy threshold of 4xl0 18 eV is, by coin¬ 
cidence, close to the value obtained for the ankle, (5.25 ± 0.12) xlO 18 eV, in the combined 
analysis of hybrid events and SD events [9]. It is thus not possible to analyse the ankle tran¬ 
sition properly. Nevertheless, the raw data do clearly display such a feature (see figure 3) 
which we do not analyse any further here due to the large uncertainty in the efficiency in 
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Figure 6. Systematic uncertainties in the derived flux as a function of the cosmic-ray energy. For a 
description of the different sources of systematic uncertainty, see the text. 


Parameter 

Power laws 

Power laws + smooth suppression 

72 {E > -Eankle) 

2.71 ± 0.02 ± 0.1 (sys) 

2.70 ±0.02 ±0.1 (sys) 

T-break 

(4.01 ± 0.21±};° (sys)) x 10 19 eV 


73 (E > T^break) 

E s 

5.98 ± 0.61 (sys) 

(5.12 ± 0.25+})] (sys))xl0 19 eV 

A 7 


5.4 ± l.Ol};} (sys) 

X 2 /ndof 

15.7/10 

13.2/10 


Table 1. Fitted parameters, with statistical and systematic uncertainties, parameterising the energy 
spectrum measured with the inclined events. 


this region (see figure 2 ). As a consequence, the 71 and Rankle parameters in equation (4.1) 
have been fixed to the values reported in [9]. 

Both models consistently describe the “flat” region of the spectrum above the ankle up 
to the observed onset of the suppression at ~4x 10 19 eV by a power law with a spectral index 
of 72 = 2.7. The spectral index after the steepening is less certain due to the low number of 
events and large systematic uncertainties. The significance of the suppression'^ is ~6.6 (t with 
220.4 events expected above the break energy inbreak while only 102 events were observed. 

A different observable that characterises the suppression is the energy £ 7/2 at which 
the integral spectrum drops by a factor of two relative to the power-law extrapolation from 
lower energies, as suggested in [38]. Here, the integral spectrum was computed by integrating 

^Following the recipe given in [37], the null hypothesis that the power law with spectral index 72 = 
2.71 ± 0.02 continues beyond the suppression point can be rejected due to the low probability (5jj®) x 10 -11 . 
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Figure 7. Energy spectrum of the cosmic rays, corrected for energy resolution, derived from inclined 
data fitted with simple power laws (dashed line) and with eqs. (4.1) and (4.2) (solid line). The 
systematic uncertainty on the energy scale is 14%. 


the parameterisation given by eqs. (4.1) and (4.2) and the parameters reported in table 1. 
The result yields £ 1/2 = (3.21 ± 0.01 ± 0.8 (sys)) x 10 19 eV, which is smaller than the value of 
~ 5.3xl0 19 eV predicted for the GZK energy cutoff of protons [38] (practically independent 
of the generation index of the assumed energy spectrum), with a difference at the level of 
2 .6(7. In [38] the assumption was made that the sources of ultra-high energy cosmic rays are 
uniformly distributed over the universe and are accelerators of protons. In reality, sources are 
discrete and in the GZK region the shape of the spectrum will be dominated by the actual 
distribution of sources around us. Further discussion of this point is given in [39]. Other 
scenarios for the high-energy suppression of the spectrum (e.g., [40-43]) attribute it to the 
limiting acceleration energy at the sources rather than to the GZK effect, providing a good 
description of the combined Auger energy spectrum as shown in [44]. 

To compare the spectrum obtained here with other measurements, we have adopted the 
technique suggested in [45] in which the differential flux at each energy is compared with 
the expected differential flux from a reference spectrum. We choose as reference a spectrum 
with an index 1 of 2.67 fitted to the flux presented here (see figure 5) in the energy bin 
corresponding to log 10 (£/eV) = 18.95 (bin width of 0.1), which contains over 1425 events. 
The reference spectrum is thus 

J ex p = 2.52xl0 31 £" 2 - 67 eV“ 1 km“ 2 sr“ 1 yr- 1 . (5.1) 

4 The index of the reference spectrum has no physical significance, so its choice is quite arbitrary. 
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Figure 8. Fractional difference between the energy spectrum of cosmic rays derived from SD data 
with 9 ^ 60° recorded at the Pierre Auger Observatory and a reference spectrum with an index of 
2.67. Residuals for spectra derived from Auger SD data with 9 < 60° [9] and the Telescope Array SD 
data with 9 < 55° [12] are also shown for comparison. 


In figure 8 the spectrum obtained with inclined events is displayed as a fractional differ¬ 
ence (also called the residual) with respect to the reference spectrum, in comparison to the 
residual for the spectrum obtained from data with zenith angles less than 60° [9] . The com¬ 
parison shows that both spectra are in agreement within errors. We note that the last two 
points are systematically below the corresponding measurement obtained for Auger events 
with zenith angles less than 60°. Also the point corresponding to log 10 (-E/eV) = 19.35 is 
below the vertical one. However, the differences are at the 2 a level. 

Figure 8 also displays the spectrum obtained from SD events with zenith angles less 
than 55° recorded with the Telescope Array (TA) detector [12]. The comparison of the 
residuals for the three spectra (also illustrated in the form of J E 3 in figure 9) shows that in 
the region between the ankle and the suppression the Auger spectra fit well due to our choice 
of reference spectrum, and the average of the residuals for TA is about +23%. The spectra 
determined by the Auger and TA Observatories are consistent in normalisation and shape 
within the systematic uncertainties in the energy scale quoted by the two collaborations [46] . 
However, differences are clearly seen in the high-energy region and are not understood thus 
far. Understanding the origin of this difference, whether from anisotropies at high energies, 
composition-related effects, experimental issues or any other reason, is of high priority in 
the efforts to understand the origin of the UHECRs. However the dedicated study of this 
discrepancy is beyond the scope of this work. Since 2012 there has been a collaborative 
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Figure 9. Compilation of the energy spectrum of cosmic rays derived from SD data recorded at the 
Pierre Auger Observatory with 9 ^ 60° (circles) and 9 < 60° (squares), and at the Telescope Array 
with 9 < 55° (triangles). 


effort between the Pierre Auger and Telescope Array Collaborations to investigate the level 
of agreement between the different measurements of the UHECR energy spectrum and to 
understand the sources of possible discrepancies, by examining the different measurement 
techniques and analysis methods employed by these groups. The latest results obtained by 
the energy spectrum working group can be found in [46]. 

6 Summary 

The cosmic-ray spectrum has been obtained for energies exceeding 4xl0 18 eV using showers 
with zenith angles between 60° and 80° recorded with the surface detector of the Pierre 
Auger Observatory during the time period between 1 January 2004 and 31 December 2013. 
It has been shown that the SD array becomes fully efficient for inclined events above this 
energy. The results can be described by a power-law spectrum with spectral index 2.7 above 
5xl0 18 eV and clearly indicate a steepening of the cosmic-ray spectrum above an energy 
around 4xl0 19 eV. These features and the normalisation of the spectrum are in agreement 
with previous measurements made with the Pierre Auger Observatory (using SD data and 
hybrid data with zenith angles less than 60°) and also with the measurements of the Telescope 
Array. 

The inclined data set is independent and complementary to the vertical data set and its 
reconstruction is performed using an independent method. These data provide a 29% increase 
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in number of events over the previous event set. In addition to obtaining an independent 
measurement of the cosmic-ray spectrum reported here, the inclined data are being analysed 
to explore primary composition, to constrain the current models that attempt to describe 
the hadronic interactions at energies above 4xl0 18 eV (i?cM ~ 87TeV), and to improve the 
studies of the arrival directions of the cosmic rays by extending the accessible fraction of the 
sky to 85%. 
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